str(cat$countries)

pays<-cat$countries[cat$countries$income_level!="Aggregates",c("iso3c", "country","capital_city","longitude","latitude", "region","income_level")]
kable(head(pays))

indic<-cat$indicators
kable(head(indic,3))

vars <- wbsearch(pattern = "GDP",fields="indicator")
kable(head(vars))

meta<-cat$indicators[cat$indicators$indicator_id %in% c("NY.GDP.MKTP.KD","NY.GDP.MKTP.PP.CD"),]
kable(meta)

tabGDP <- wb_data(indicator = c("NY.GDP.MKTP.KD"),
                 return_wide = TRUE,
                 mrv=1,
                 country ="countries_only")



tabGDP<-tabGDP[order(tabGDP$NY.GDP.MKTP.KD,decreasing = T),]
kable(head(tabGDP,10))

tabGDP <- wb_data(indicator = c("NY.GDP.MKTP.PP.CD"),
                 return_wide = TRUE,
                 mrv=1,
                 country ="countries_only")



tabGDP<-tabGDP[order(tabGDP$NY.GDP.MKTP.PP.CD,decreasing = T),]
kable(head(tabGDP,10))

world_data <- wb_data(indicator = c("AG.SRF.TOTL.K2","AG.LND.ARBL.HA","SP.POP.TOTL","SP.URB.TOTL","NY.GDP.MKTP.CD", "NY.GDP.MKTP.PP.CD"),
                 return_wide = TRUE,
                 start_date = 1960,
                 end_date = 2024,
                 country ="countries_only")

world_data<-world_data[,-c(1,3)]

# recodage (attention : ordre alphabetique)
names(world_data)<-c("iso3c","date","ARB","SRF","GDP","GPP","POP","URB")

world_data$SRF[world_data$iso3c=="CAN"]<-9984670

pays<-cat$countries[cat$countries$income_level!="Aggregates",c("iso3c", "country","capital_city","longitude","latitude", "region","income_level")]
tab<-merge(pays,world_data, by="iso3c",all.x=F,all.y=T)
kable(head(tab))

saveRDS(tab,"data/wb_don_1960_2023.Rdata")

# Extract meta
meta<-cat$indicators[cat$indicators$indicator_id %in% c("AG.SRF.TOTL.K2","AG.LND.ARBL.HA","SP.POP.TOTL","SP.URB.TOTL","NY.GDP.MKTP.CD", "NY.GDP.MKTP.PP.CD"),]

# Select information and add personal code
meta<-data.frame(meta[,c(1,2,4,5)])
meta$shortcode<-c("ARB","SRF","GDP","GPP","POP","URB")
meta<-meta[,c(5,1,2,3,4)]

# Display
kable(meta)

# Save meta
saveRDS(meta,"data/wb_don_1960_2023_meta.Rdata")
write.table(meta,"data/wb_don_1960_2023_meta.csv",sep=";",dec=",", row.names = F,fileEncoding = "UTF-8")

# Load data
don<-readRDS("data/wb_don_1960_2023.Rdata")

# Correction terres arables
length(don$ARB[don$date==2016])

don$ARB[don$date==2022]<-don$ARB[don$date==2021]
don$ARB[don$date==2023]<-don$ARB[don$date==2021]
don$ARB[don$date==1960]<-don$ARB[don$date==1961]


## Correction superfice
don$SRF[don$date==2023]<-don$SRF[don$date==2021]
don$SRF[don$date==1960]<-don$SRF[don$date==1961]

tabmis<-don %>% select(iso3c,country,date,ARB,SRF,GDP,GPP,POP, URB)  %>%
  gather("ARB", "SRF","GDP","GPP","POP","URB",key = VAR, value = DON) %>%
  mutate(Missing = is.na(DON))

kable(head(tabmis))             

res <- tabmis %>% group_by(VAR, date) %>% 
                  summarise(nbmis=sum(Missing), nb=n(), pct = 100*nbmis/nb) 

p<-ggplot(res) + 
  aes(x=date,y=pct,color=VAR) +
  geom_line() +
  scale_y_continuous("% de pays à valeurs manquantes", breaks = c(0,10,20,30,40,50,60,70,80,90,100))+
  scale_x_continuous("Année", breaks = c(1960, 1965, 1970, 1975, 1980, 1985,1990, 1995, 2000, 2005, 2010, 2015, 2020, 2025))
p

res <- tabmis %>% group_by(iso3c, country) %>% 
                  summarise(nbmis=sum(Missing), nb=n(), pct = 100*nbmis/nb) %>%
                  arrange(-pct)

kable(head(res,20))

# load a data.frame with columns space, time, indic 
estim<-function(df   = don)               # dataframe
 
{


# Select variable
sel<-dcast(df,space~time)
M<-as.matrix(sel[,-1])
rownames(M)<-sel$space

####################################################
########   Estimation Procedure ####################
####################################################

# check for lines or column  filled with missing values and eliminate
check1<-apply(!is.na(M),FUN="sum",1)
check2<-apply(!is.na(M),FUN="sum",2)
M<-M[check1>0,check2>0]



dim(M)


# define dimensions
nr<-nrow(M)
nr
nc<-ncol(M)
nc


# create time matrix with NA
Mt<-matrix(rep(1:nc,nr),nrow=nr,nc=nc,byrow=T)
rownames(Mt)<-rownames(M)
colnames(Mt)<-colnames(M)
Mt[is.na(M)]<-NA
Mt[1:4,1:4]



# identify previous and next available valuer
Mt_min<-Mt
Mt_max<-Mt
for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if (is.na(Mt[i,j])) {
      Mt_min[i,j]<-max(Mt[i,1:j],na.rm=T)
      Mt_max[i,j]<-min(Mt[i,j:nc],na.rm=T)    
    }  
  }
}
Mt_min[1:4,1:4]
Mt_max[1:4,1:4]


# choose estimation method
M_met<-matrix("OK",nrow=nr,ncol=nc)
rownames(M_met)<-rownames(M)
colnames(M_met)<-colnames(M)
for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if (is.na(Mt[i,j])) {
      M_met[i,j]<-"IN"
      if (is.infinite(Mt_max[i,j])) {M_met[i,j]<-"EX"} 
      if (is.infinite(Mt_min[i,j])) {M_met[i,j]<-"RE"}  
    }   
  }
}
M_met[1:4,1:4]

###### Estimation of missing values  ####

M_est<-M
str(M_est)

#### step 1: Interpolation ######

for ( i in 1:nr) { 
  for ( j in 1:nc) { 
    if ((M_met[i,j]=="IN")) {
      t0<-Mt_min[i,j]
      t1<-Mt_max[i,j] 
      tacm<-(M[i,t1]/M[i,t0])**(1/(t1-t0))
      M_est[i,j]<-M[i,t0]*(tacm**(j-t0))
    }   
  }
}
M_est[1:4,1:4]


#### step 2: Extrapolation ######

for ( i in 1:nr) { 
  for ( j in 2:nc) { 
    if ((M_met[i,j]=="EX")) {
      N<-M_est
      t0<-j-1
      t1<-j
      N<-M_est[is.na(M_est[,t0])==F,]
      N<-N[is.na(N[,t1])==F,]
      tacm<-(sum(N[,t1])/sum(N[,t0]))**(1/(t1-t0)) 
      M_est[i,j]<-M_est[i,t0]*(tacm**(t1-t0))
    }   
  }
}
M_est[1:4,1:4]


#### step 3: Retropolation ######

for ( i in 1:nr) { 
  for ( j in (nc-1):1) { 
    if ((M_met[i,j]=="RE")) {
      N<-M_est
      t0<-j
      t1<-j+1
      N<-M_est[is.na(M_est[,t0])==F,]
      N<-N[is.na(N[,t1])==F,]
      tacm<-(sum(N[,t1])/sum(N[,t0]))**(1/(t1-t0)) 
      M_est[i,j]<-M_est[i,t1]/(tacm**(t1-t0))
    }   
  }
}
M_est[1:4,1:4]



#######################################################
############# EXPORT RESULTS  #########################
#######################################################

x<-reshape2::melt(M_est)
names(x)<-c("space","time","estim")
y<-reshape2::melt(M_met)
names(y)<-c("space","time","method")
z<-merge(x,y, by=c("space","time"))

return(z)
}

# Estim ARB
df<-data.table(space=don$iso3c,time=don$date,index=don$ARB)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","ARB_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim SRF
df<-data.table(space=don$iso3c,time=don$date,index=don$SRF)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","SRF_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)


# Estim GDP
df<-data.table(space=don$iso3c,time=don$date,index=don$GDP)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","GDP_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

## Estim GPP
df<-data.table(space=don$iso3c,time=don$date,index=don$GPP)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","GPP_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)


# Estim POP
df<-data.table(space=don$iso3c,time=don$date,index=don$POP)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","POP_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

# Estim URB
df<-data.table(space=don$iso3c,time=don$date,index=don$URB)
est<-estim(df)[,1:3]

names(est)<-c("iso3c","date","URB_est")
don<-merge(don,est, by=c("iso3c","date"), all.x=T,al.y=F)

obs<-don %>% filter(is.na(GPP_est)==F) %>%select(iso3c, date, POP_est, GDP_est, GPP_est)


obs$y<-obs$GPP_est/obs$POP_est
obs$x<-obs$GDP_est/obs$POP_est
plot(obs$x,obs$y,log = 'xy')

lm(log(obs$y)~log(obs$x))

mod<-lm(log(obs$GPP_est)~log(obs$GDP_est)+log(obs$POP_est))
summary(mod)

# Estimation
don$GPP_est[is.na(don$GPP_est)]<-don$POP_est[is.na(don$GPP_est)]*((don$GDP_est[is.na(don$GPP_est)]/don$POP_est[is.na(don$GPP_est)])**0.758)*exp(2.745)



# Save
saveRDS(don,"data/wb_don_1960_2023.Rdata")

don<-readRDS("data/wb_don_1960_2023.Rdata")
#don<-as_tibble(don)
don2<-don %>%  group_by(date) %>% mutate(ARB_rnk = rank(desc(ARB_est)),
                                         SRF_rnk = rank(desc(SRF_est)),
                                         GDP_rnk = rank(desc(GDP_est)),
                                         GPP_rnk = rank(desc(GPP_est)),                                     
                                         POP_rnk = rank(desc(POP_est)),
                                         URB_rnk = rank(desc(URB_est)),
                                         TOT_rnk = (ARB_rnk+SRF_rnk+GDP_rnk+GPP_rnk+POP_rnk+URB_rnk)/6) 

tab<-don2 %>% filter(date==1960) %>% 
             select(iso3c, country, TOT_rnk, ARB_rnk,SRF_rnk,GDP_rnk,GPP_rnk,POP_rnk,URB_rnk,TOT_rnk) %>%
             arrange(TOT_rnk)
kable(head(tab,20), digits=1)

tab<-don2 %>% filter(date==2023) %>% 
             select(iso3c, country, TOT_rnk, ARB_rnk,SRF_rnk,GDP_rnk,GPP_rnk,POP_rnk,URB_rnk,TOT_rnk) %>%
             arrange(TOT_rnk)
kable(head(tab,20), digits=1)

don2<-don2 %>%  group_by(date) %>% mutate(ARB_wgt = 100*ARB_est/sum(ARB_est,na.rm=T),
                                         SRF_wgt = 100*SRF_est/sum(SRF_est,na.rm=T),
                                         GDP_wgt = 100*GDP_est/sum(GDP_est,na.rm=T),
                                         GPP_wgt = 100*GPP_est/sum(GPP_est,na.rm=T),
                                         POP_wgt = 100*POP_est/sum(POP_est,na.rm=T),
                                         URB_wgt = 100*URB_est/sum(URB_est,na.rm=T),
                                         TOT_wgt = (ARB_wgt+SRF_wgt+GDP_wgt+GPP_wgt+POP_wgt+URB_wgt)/6)

tab<-don2 %>% filter(date==1960) %>% 
             select(iso3c, country, TOT_wgt, ARB_wgt,SRF_wgt,GDP_wgt,GPP_wgt,POP_wgt,URB_wgt,TOT_wgt) %>%
             arrange(desc(TOT_wgt))
kable(head(tab,20), digits=1)

tab<-don2 %>% filter(date==2023) %>% 
             select(iso3c, country, TOT_wgt, ARB_wgt,SRF_wgt,GDP_wgt,GPP_wgt,POP_wgt,URB_wgt,TOT_wgt) %>%
             arrange(desc(TOT_wgt))
kable(head(tab,20), digits=1)

saveRDS(don2,"data/wb_don_1960_2023.Rdata")
write.table(don2,"data/wb_don_1960_2023.csv",sep=";",dec=",", row.names = F,fileEncoding = "UTF-8")

map<-st_as_sf(countries110)
par(mar=c(0,0,0,0), mfrow=c(1,1))
plot(map$geometry,col="lightyellow")

dim(map)

sel<-map[map$adm0_a3 %in% c("FRA", "NCL","ATA","ATF","USA", "PRI","CHN","TWN","MAR", "SAH","CHN","TWN","ISR","PSX"),c("sovereignt","sov_a3","type","admin", "adm0_a3","name","note_adm0","iso_a3","wb_a3")]
kable(sel)

map<-st_as_sf(sovereignty110)
par(mar=c(0,0,0,0))
plot(map$geometry,col="lightyellow")

dim(map)

par(mfrow=c(1,2))
plot(map[map$sov_a3 %in% c("IS1","JOR","SYR","LBN","EGY"),]$geometry, col=c("gray80","orange","gray80","gray80","gray80"))
title("Limits of Israël",cex=0.5)
plot(map[map$sov_a3 %in% c("MAR","SAH","DZA","MRT"),]$geometry, col=c("gray70","orange","gray70","lightyellow"))
title("Limits of Morocco")

map<-st_as_sf(countries110)
small<-st_as_sf(tiny_countries110)
par(mar=c(0,0,0,0), mfrow=c(1,1))
plot(map$geometry,col="lightyellow")
plot(small$geometry,col="red", add=T)

map<-st_as_sf(countries50)

par(mar=c(0,0,0,0))
plot(map$geometry,col="lightyellow")

# Load map
map<-st_as_sf(countries110)
map<-st_transform(map, "+proj=laea +x_0=0 +y_0=0 +lon_0=0 + lat_0=40")
plot(map$geometry)

map<-map[c("adm0_a3","name")]
names(map)<-c("iso3c","name","geometry")

# Add polygons center
coo<-st_coordinates(st_centroid(map,of_largest_polygon = T))
map$Lon<-coo[,1]
map$Lat<-coo[,2]


# adjust some codes
map$iso3c[map$iso3c=="KOS"]<-"XKX"    # Kosovo
map$iso3c[map$iso3c=="PSX"]<-"PSE"    # Palestinian territories
map$iso3c[map$iso3c=="SDS"]<-"SSD"    # South Sudan

# Save
st_write(map,"data/wb_map_low.shp",delete_dsn=T)

saveRDS(map,"data/wb_map_low.Rdata")

# Load map
map<-st_as_sf(countries50)
map<-map[c("adm0_a3","name")]
names(map)<-c("iso3c","name","geometry")

# Add polygons center
coo<-st_coordinates(st_centroid(map,of_largest_polygon = T))
map$Lon<-coo[,1]
map$Lat<-coo[,2]


# adjust some codes
map$iso3c[map$iso3c=="KOS"]<-"XKX"    # Kosovo
map$iso3c[map$iso3c=="PSX"]<-"PSE"    # Palestinian territories
map$iso3c[map$iso3c=="SDS"]<-"SSD"    # South Sudan

# Save
st_write(map,"data/wb_map_high.shp",delete_dsn=T)

saveRDS(map,"data/wb_map_high.Rdata")

don<-readRDS("data/wb_don_1960_2023.Rdata")
tab<-don %>% filter(date > 1960, date < 2023) %>%
  select (iso3c, country, date, GDP_est,GDP_rnk, GDP_wgt) %>% 
  filter(iso3c %in% c("USA","CHN","JPN","RUS","IND","DEU","BRA","FRA","GBR","SAU")) 

tab<-as.data.frame(tab)


p <- plot_ly(tab,
             x = ~date,
             y = ~GDP_wgt,
             color= ~country) %>%
             add_lines()%>%
         layout(title = "PIB en dollars (1960-2023)",
         yaxis = list(title = "% de l richesse mondiale", type="log"), 
         xaxis = list(title = "Année"))
p

don<-readRDS("data/wb_don_1960_2023.Rdata")
tab<-don %>% filter(date > 1960, date < 2023) %>%
  select (iso3c, country, date, GPP_est,GPP_rnk, GPP_wgt) %>% 
  filter(iso3c %in% c("USA","CHN","JPN","RUS","IND","DEU","BRA","FRA","GBR","SAU")) 

tab<-as.data.frame(tab)


p <- plot_ly(tab,
             x = ~date,
             y = ~GPP_wgt,
             color= ~country) %>%
             add_lines()%>%
         layout(title = "PIB en dollars (1960-2023)",
         yaxis = list(title = "% de l richesse mondiale", type="log"), 
         xaxis = list(title = "Année"))
p

# Load map
map<-readRDS("data/wb_map_high.Rdata")

# Load don for 2015
don<-readRDS("data/wb_don_1960_2023.Rdata")
don<-don[don$date==2023,]

# Merge with map
mapdon<-merge(don,map,by="iso3c",all.x=F,all.y=T)

# Create map



g <- list(showframe = TRUE,
          framecolor= toRGB("gray50"),
          coastlinecolor = toRGB("black"),
          showland = TRUE,
          landcolor = toRGB("lightyellow"),
          showcountries = TRUE,
          countrycolor = toRGB("black"),
          countrywidth = 0.2,
         projection = list(type = 'azimuthal equal area')
 #         projection = list(type = 'Mercator')
          )



p<- plot_geo(mapdon)%>%
  add_markers(x = ~Lon,
              y = ~Lat,
              sizes = c(0, 1000),
              size = ~GDP_wgt,
              color= "red",
              hoverinfo = "text",
              text = ~paste('Pays: ', country,
                            '<br /> % mondial : ',round(GDP_wgt,3),
                            '<br /> rang mondial : ',round(GDP_rnk,0))) %>%
  layout(geo = g,
         title = "Part du PIB mondial (ppa) en 2023",
         width=800,
         height = 400)
p